#Clear R's memory as a good habit
rm(list=ls())

#Specify your working directory in R
setwd("C:/dataR/R2022ConsCourt")

library(haven)
library(survival)

mydata <- read_dta("KimNolettePOPAnalysisStata12.dta")
summary(mydata)

#Note that our R codes for Models 4 and 5 in Table 1 below rely on Box-Steffensmeier, De Boef, and Joyce's (2007).
##############################################################################
##Robustness Check: Model 4 in Table 1: Conditional Frailty Model           ##
##############################################################################
condfrailty.gamma.em <- coxph(Surv(istart,ifinish,event)~
                  commonlaw+VCengage+
                  polinsur+RoverG+dt3yrs+v2xnp_pres12+lnlatent_gdppcWB22+
                  strata(eventnum)+frailty(id,distribution="gamma",method="em"),
                  method=c("efron"),data=mydata)
summary(condfrailty.gamma.em)

nofrailty.gamma.em <- coxph(Surv(istart,ifinish,event)~
                  commonlaw+VCengage+
                  polinsur+RoverG+dt3yrs+v2xnp_pres12+lnlatent_gdppcWB22+
                  strata(eventnum),method=c("efron"),data=mydata)
summary(nofrailty.gamma.em)

#################################
#Obtain Log-Likelihood for Model#
#################################
condfrailty.gamma.em.loglike2 <- condfrailty.gamma.em$loglik[[2]]
condfrailty.gamma.em.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
condfrailty.gamma.em.theta <- condfrailty.gamma.em$history[[1]]$theta
condfrailty.gamma.em.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
condfrailty.gamma.em.ill <- condfrailty.gamma.em$history[[1]]$c.loglik
condfrailty.gamma.em.ill

#Obtain likelihood for the no frailty model
nofrailty.loglike2 <- nofrailty.gamma.em$loglik[[2]]
nofrailty.loglike2

#Calculate LR test for theta, report this test for whether theta is significant
#nofrailty.gamma.em.loglike2 - condfrailty.gamma.em.ill
#The numbers in parentheses below come from the two estimates above: 307.3081, 307.3081
#Remove minus(-) from both before subtraction
condfrailty.gamma.em.lr <- 2*(307.3081-307.3081)
condfrailty.gamma.em.lr

#Calculate p-value, it has 1 df
condfrailty.gamma.em.lrtest <- 1-pchisq(condfrailty.gamma.em.lr,df=1)
condfrailty.gamma.em.lrtest


##############################################################################
##Robustness Check: Model 5 in Table 1: Alternative Dependent Variable      ##
##############################################################################
condfrailty.alt.DV <- coxph(Surv(istart_ccpcnc,ifinish_ccpcnc,event_ccpcncyr1)~
                  commonlaw+VCengage+
                  polinsur+RoverG+dt3yrs+v2xnp_pres12+lnlatent_gdppcWB22+
                  strata(eventnum_ccpcnc)+frailty(id,distribution="gamma",method="em"),
                  method=c("efron"),data=mydata)
summary(condfrailty.alt.DV)

nofrailty.alt.DV <- coxph(Surv(istart_ccpcnc,ifinish_ccpcnc,event_ccpcncyr1)~
                  commonlaw+VCengage+
                  polinsur+RoverG+dt3yrs+v2xnp_pres12+lnlatent_gdppcWB22+
                  strata(eventnum_ccpcnc),method=c("efron"),data=mydata)
summary(nofrailty.alt.DV)

#################################
#Obtain Log-Likelihood for Model#
#################################
condfrailty.alt.DV.loglike2 <- condfrailty.alt.DV$loglik[[2]]
condfrailty.alt.DV.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
condfrailty.alt.DV.theta <- condfrailty.alt.DV$history[[1]]$theta
condfrailty.alt.DV.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
condfrailty.alt.DV.ill <- condfrailty.alt.DV$history[[1]]$c.loglik
condfrailty.alt.DV.ill

#Obtain likelihood for the no frailty model
nofrailty.loglike2 <- nofrailty.alt.DV$loglik[[2]]
nofrailty.loglike2

#Calculate LR test for theta, report this test for whether theta is significant
#nofrailty.alt.DV.loglike2 - condfrailty.alt.DV.ill
#The numbers in parentheses below come from the two estimates above: 363.5052, 363.5052
#Remove minus(-) from both before subtraction
condfrailty.alt.DV.lr <- 2*(363.5052-363.5052)
condfrailty.alt.DV.lr

#Calculate p-value, it has 1 df
condfrailty.alt.DV.lrtest <- 1-pchisq(condfrailty.alt.DV.lr,df=1)
condfrailty.alt.DV.lrtest

##############################################################################
##Robustness Check: Model 1 in Table OA3: Fine & Gray Model                 ##
##############################################################################
library(cmprsk)
library(crrSC)
library(aod)

cov1=cbind(mydata$commonlaw,mydata$VCengage,mydata$polinsur,mydata$RoverG,
           mydata$dt3yrs,mydata$v2xnp_pres12,mydata$lnlatent_gdppcWB22)
finegray <- crr(ftime=mydata$time_CR_ccpcnc, fstatus=mydata$event_CR_ccpcncyr1, 
            cov1, failcode=2, cencode=0, 
            na.action=na.omit, gtol=1e-6, maxiter=10, init=0)
summary(finegray)
finegray$loglik
wald.test(b=finegray$coef, Sigma=finegray$var, Terms=1:7)

##############################################################################
##Robustness Check: Model 2 in Table OA3: Stratified Fine & Gray Model      ##
##############################################################################

cov1=cbind(mydata$commonlaw,mydata$VCengage,mydata$polinsur,mydata$RoverG,
           mydata$dt3yrs,mydata$v2xnp_pres12,mydata$lnlatent_gdppcWB22)
finegraystrata <- crrs(ftime=mydata$time_CR_ccpcnc, fstatus=mydata$event_CR_ccpcncyr1, 
                  cov1, strata=mydata$eventnum_ccpcnc, failcode=2, cencode=0,
                  ctype=1, na.action=na.omit, gtol=1e-6, maxiter=10, init=0)
print(finegraystrata)
finegraystrata$loglik
wald.test(b = finegraystrata$coef, Sigma = finegraystrata$var, Terms = 1:7)

